use rand::Rng;

// 求最大公约数
fn gcd(a: u128, b: u128) -> u128 {
    if b == 0 {
        a
    } else {
        gcd(b, a % b)
    }
}

// 快速乘 (带模乘法)
fn bmul(a: u128, b: u128, m: u128) -> u128 {
    let a = a % m;  
    let b = b % m;  
    let c = a.checked_mul(b).unwrap_or(m);  
    let result = c % m;  
    result
}


/// 快速幂计算 a^b % m
pub fn quick_pow(mut a: u128, mut b: u128, m: u128) -> u128 {
    let mut result = 1;
    a %= m;
    while b > 0 {
        if b & 1 == 1 {
            result = result * a % m;
        }
        a = bmul(a, a, m);
        b >>= 1;
    }
    result
}
/// 米勒-拉宾素性测试
fn miller_rabin(n: u128, test_times: usize) -> bool {
    if n < 2 { return false; }
    if n == 2 || n == 3 { return true; }
    if n % 2 == 0 || n % 3 == 0 { return false; }

    let mut u = n - 1;
    let mut t = 0;
    while u % 2 == 0 {
        u >>= 1;
        t += 1;
    }

    let mut rng = rand::thread_rng();
    for _ in 0..test_times {
        let a = rng.gen_range(2..=n-2);
        let mut v = quick_pow(a, u, n);
        if v == 1 || v == n - 1 { continue; }

        let mut s = 0;
        while s < t - 1 {
            v = quick_pow(v, 2, n);
            if v == n - 1 { break; }
            s += 1;
        }
        if v != n - 1 { return false; }
    }
    true
}

/// 综合素性测试函数
pub fn is_prime(n: u128) -> bool {
    if n <= 1 {
        return false;
    }
    if n <= 3 {
        return true;
    }
    if n % 2 == 0 || n % 3 == 0 {
        return false;
    }
    if n < 1000000 {
        // 对于较小的数，使用试除法
        let mut i: u128 = 5;
        while i * i <= n {
            if n % i == 0 || n % (i + 2) == 0 {
                return false;
            }
            i += 6;
        }
        return true;
    } else {
        // 对于较大的数，使用米勒-拉宾素性测试
        return miller_rabin(n, 8); // 测试次数可以根据需要调整
    }
}

// Pollard-Rho 算法
fn pollard_rho(x: u128) -> u128 {
    let mut rng = rand::thread_rng();
    let mut s = 0;
    let mut t = 0;
    let c = rng.gen_range(1..x);
    let mut val = 1;

    // 倍增优化
    for goal in (1..).map(|goal| goal * 2) {
        for step in 1..=goal {
            t = (bmul(t, t, x) + c) % x;
            val = bmul(val, (t as i128 - s as i128).abs() as u128, x);
            if step % 127 == 0 {
                let d = gcd(val, x);
                if d > 1 {
                    return d;
                }
            }
        }

        let d = gcd(val, x);
        if d > 1 {
            return d;
        }

        s = t;
        val = 1;
    }

    x
}

// 分解数
fn fac(mut x: u128, max_factor: &mut u128) {
    if x <= *max_factor || x < 2 {
        return;
    }

    // 如果是质数，直接更新最大质因子
    if is_prime(x) {
        *max_factor = (*max_factor).max(x);
        return;
    }
    let mut  p = x;
    while p >= x {
        p = pollard_rho(x);
    }

    // 分解掉因子 p
    while x % p == 0 {
        x /= p;
    }

    // 递归分解剩余部分
    fac(x, max_factor);
    fac(p, max_factor);
}

// 查找最大质因子
pub fn find_max_prime_factor(number: u128) -> u128 {
    let mut max_factor = 0;
    fac(number, &mut max_factor);
    if max_factor == number {
        number
    } else {
        max_factor
    }
}